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We describe a method for simulating the real time evolution of extended quantum systems in two 
dimensions. The method combines the benefits of integrability and matrix product states in one 
dimension to avoid several issues that hinder other applications of tensor based methods in 2D. In 
particular it can be extended to infinitely long cylinders. As an example application we present 
results for quantum quenches in the 2D quantum (2+1 dimensional) Ising model. In quenches that 
cross a phase boundary we find that the return probability shows non-analyticities in time. 


The advent of ultra cold atomic gas experiments has 
led to a surge of interest in the time evolution and out-of- 
equilibrium behaviour of many-body quantum systems. 
Much effort has been focused on one dimensional (ID) 
problems because these can be tackled by analytically 
tractable or highly accurate numerical methods. Key 
questions that these studies have sought to elucidate are 
whether and how such systems thermalise after a sudden 
change, or ‘quantum quench’ of a system’s Hamiltonian; 
with particular emphasis on the role played by conserved 
charges in ID integrable systems [1 11]. 

Experiments however, are not limited to ID and it is 
interesting to explore similar questions in two dimensions 
(2D) and above [12]. Unfortunately there is no analogue 
in 2D of the aforementioned analytically exact ID meth¬ 
ods. Numerical approaches using matrix product state 
(MPS) representations, so successful in ID, suffer in 2D 
due to the ‘area law’ growth of entanglement [13, 14], 
This growth reduces the efficiency of MPS (and related 
‘tensor’) algorithms and limits them to smaller system 
sizes. 

Nonetheless MPS algorithms can be applied in 2D, by 
labeling lattice sites (usually in a zigzag fashion) to map 
to a ID system [15]. The cost is that nearest neighbor in¬ 
teractions in 2D are mapped to increasingly long ranged 
ID interactions, imposing an increasing numerical bur¬ 
den. Recently progress has been made in performing real 
time evolution on MPS with such long ranged Hamilto¬ 
nians by two different routes [16, 17]. Algorithms based 
on generalizations of MPS to higher dimensions, such as 
projected entangled pair states (PEPS) [18, 19], make 
use of imaginary time evolution to find ground states 
[20]. However these higher dimensional tensor methods 
have not been applied to real time evolution. 

In this letter we demonstrate that real time evolution 
is possible for large 2D systems by combining informa¬ 
tion coming from exactly solvable models with a highly 
anisotropic MPS formulation. Such an approach retains 
the contraction efficiency of matrix product states over 
other tensor methods, while avoiding the build up of 



FIG. 1. Anisotropic setup for a 2D system as an array of 
N chains of length R, coupled by an interaction J±. The 
cylinder can be joined together at its ends to study toroidal 
systems. 

long ranged interactions. Our setup will be similar to 
that used in the density matrix renormalisation group 
(DMRG) studies described in Refs. [21, 22] except that 
here we are explicit in our use of MPS. This change allows 
for straightforward implementation of algorithms other 
than DMRG, including those for time evolution and for 
accurately working with the thermodynamic limit. In 
particular using time evolving block decimation (TEBD) 
[23] we demonstrate that we can study the time evolu¬ 
tion after a quench of infinitely long cylinders, with suf¬ 
ficient circumference that we approach the 2D thermo¬ 
dynamic limit. This includes strong quenches where we 
cross phase boundaries of a 2D quantum system. 
Method: At the core of our method is the wish to max¬ 
imise the analytically exact input going into our MPS 
algorithm, while simultaneously controlling the growth 
of entanglement entropy. The construction we use is de¬ 
picted in Fig. 1: a coupled array of exactly solvable ID 
subunits. For each subunit, we have exact knowledge of 
the spectrum and matrix elements. This exact knowledge 
means that we begin with the numerics already having 
accounted for much of the strong correlations of the sys¬ 
tem. We emphasize our use of exactly solvable models 
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as a building block is not much of a limitation to the 
method. Such models are ubiquitous in ID, including 
Heisenberg spin chains, Luttinger liquids, and Hubbard 
models to name but a few [24, 25]. In this framework, a 
state of a system of N chains is written in MPS form via 

\*) 2D = 'ZA°M--A™W\ <Tl ... aN ) (!) 


where each matrix A ai M is labelled by a chain i and an 
eigenstate of that individual chain ct,; . Like the single sites 
used in ID MPS algorithms, we are able to manipulate 
these chain eigenstates because we know their energies 
and matrix elements for any relevant operator. 

For ground and low-lying states of the system the en¬ 
tanglement entropy Se scales as the boundary ‘area’, 
that is to say the chain length. By keeping the chain 
length finite we can throttle the growth of Se- By part¬ 
nering this with the fact that for the systems that we 
will study, finite size effects are exponentially suppressed , 
we are able to keep Se small while remaining in the 2D 
thermodynamic limit. We have previously demonstrated 
the effectiveness of this methodology in equilibrium by 
studying a 2D quantum (i.e. 2 + 1 dimensional) critical 
point [21, 22]. 

The continuum ID subunits will necessarily have an 
infinitely large Hilbert space. However if the system size 
R is finite the spectrum is discrete, and we may trun¬ 
cate at a cutoff energy E c . This step is justified by ap¬ 
peal to the truncated conformal spectrum approach [26] 
where it has been observed over a wide body of examples 
[27-29] that for relevant (in the renormalisation group 
sense) interchain interactions, the low energy sector of 
a perturbed integrable system is formed primarily from 
(possibly strong) admixtures of low lying states of the 
unperturbed system. Here we will focus on exactly such 
interchain perturbations. 

Eq. 1 differs from a MPS for a ID system only in 
that the ‘physical indices’ a may be large (see Table I in 
[30]), requiring strict use of sparse matrices to maximise 
computational resources. It is also important to take ad¬ 
vantage of good quantum numbers and to perform ma¬ 
trix operations (e.g. singular value decompositions) in a 
block diagonal manner, to help preserve the sparse nature 
of the matrices and increase numerical efficiency. 

MPS time evolution algorithms may then be imple¬ 
mented just as for a ID system, including TEBD [23] 
and its infinite counterpart (iTEBD) [31, 32], For the 
former we may work with a torus or open cylinder ge¬ 
ometry; the latter corresponds to an infinitely long cylin¬ 
der. Both algorithms decompose the time evolution oper¬ 
ator exp [—iHt\ into a product of N t time step operators, 
t, = N t T. Each step is itself approximately decomposed 
into a product of two site (or chain) operations. The er¬ 
ror at each step is proportional to the time increment r 
raised to a power given by the order of the decomposition. 


A more important source of error is the compression of 
the MPS after each step via Schmidt decompositions. We 
compress by fixing a minimum singular value size, s m i n : 
singular values smaller than this threshold value are dis¬ 
carded. In this sense our algorithm is adaptive, as y, 
and the degree of encoded entanglement can grow. ‘Lieb- 
Robinson’ type arguments limit the rate of growth of Se 
after a quench [33-35], but y may grow exponentially, 
limiting the maximum timescales that can be reached. 

For our 2D algorithm, forming the time evolution op¬ 
erator requires the exponentiation of a two chain Hamil¬ 
tonian, which in turn necessitates the diagonalisation of 
the same object. This is a numerically costly step, but 
need only be done once at the beginning, and the result 
stored for later use. 

In this letter we present results for quenches in the 2D 
quantum Ising model: 


#2DQI = E 


f R 

-ffiD.i + J± dx a-(x)a- +1 (x) 
Jo 


( 2 ) 


We represent the model as ID Ising chains (of index i and 
length R ) coupled together with a longitudinal spin-spin 
interaction. We take each chain ffio,i to be the con¬ 
tinuum limit of the ID lattice quantum Ising model—or 
transverse field Ising model (TFIM)- with Hamiltonian, 
- J \\ E;K/ cr ++i + (1 + g)(xfi)] with l an index along 
the chain. In the continuum limit this reduces to a the¬ 
ory of a ID Majorana field with mass A = gj ||. An¬ 
alytic expressions for the spectrum of this theory and 
the spin matrix elements are detailed in Ref. [36]; we 
summarize the salient features in [30]. Expanding the 
Majorana held in terms of fermionic modes ipt and 
ipki (the continuum versions of the usual Jordan-Wigner 
lattice fermions) yields a quadratic chain Hamiltonian 
-ffiD.i = E/w with dispersion e ki = a/A 2 + kf. 

We work in units such that the intrachain velocity, v, is 
dimensionless and equal to unity. We also define a dimen¬ 
sionless interchain coupling j_l = | A | 1 ^. For disor¬ 
dered (A < 0) chains a finite value of the interchain cou¬ 
pling jx leads to a 2D quantum (d=2+l) order-disorder 
transition at a critical value j± = j c = 0.185 [22]. 

We compute the evolution of the postquench state us¬ 
ing iTEBD and TEBD, with first and second order Trot¬ 
ter decompositions of the time evolution operator, and 
time steps r. The error associated with such decomposi¬ 
tions is dependent on j± and r, but even for the strongest 
quenches presented in this work we can choose r small 
enough for convergence (see the supplementary material 
[30]). For each set of parameters, we first establish that 
the numerical results are converged in s m i n or y before 
increasing the cutoff E c . Convergence of the method in 
Smin is demonstrated in [30]. We have also checked the 
algorithm for two analytically tractable cases: the per¬ 
turbative limit (j± <C 1) and a model of free fermionic 
chains with interchain hopping. In both cases we find 



3 



3 4 


FIG. 2. Fermion occupation number, rii(x) scaled by inter¬ 
chain coupling, jj_. We indicate the time scale tn at which 
we expect the system postquench to see the effects of the fi¬ 
nite circumference of the system. Inset: R = 10 iTEBD data 
compared with the perturbative result (P.T.) (dashed line). 



FIG. 3. Fermion occupation number, rii(x), scaled by inter¬ 
chain coupling, j± = 0.2, squared. Curves for different E c 
are shown, corresponding to more than doubling the num¬ 
ber of retained states in the chain spectrum. The agreement 
is excellent until the latest times, even though this quench 
crosses a critical point. Inset: the nearest neighbor spin-spin 
correlation function showing scaling with j± and R. 

excellent agreement with our numerical results [30]. 
Results: In the following we present results of quantum 
quenches where the initial state of system corresponds to 
the j± = 0 ground state, whereupon at t = 0 we turn 
on a finite interchain coupling j±. We focus mainly on 
results for infinitely long cylinders, leaving a discussion 
of the effect of finite chain number, N, until the end. We 
first address the question of what time scales we expect 
to feature in the quench. To provide a partial answer 


we turn to the quasiparticle causality picture of Refs. 
[1, 2, 33]. The energy imparted by the quench produces 
quasiparticle excitations which are entangled on a length 
scale |A| —1 along the chain. Intrachain scattering then 
only has an effect after a time, £a = (2i> |A|) -1 . On 
the other hand, the time scale governing interchain scat¬ 
tering can be estimated using Fermi’s golden rule to be 
tj ± = |A|( Jj_R)~ 2 . The final time scale of import 
is that encoding the chain length, R. This scale, given 
by tR ~ R/2v = |A|i?tA, describes the time for two 
quasiparticles, created at the same point and moving in 
opposite directions, to travel around a chain and then 
meet again. Hence there is a region, < t < £r, 

where we may expect the time evolution to be represen¬ 
tative of the 2D thermodynamic limit. But for t > tn 
the finite nature of the chains’ circumferences will play 
a role. We stress that £r does not govern the time scale 
for revivals in the system. Instead these occur on a much 
longer time scale, t rev i va i ~ Ntj ± where N is the number 
of chains in the system. Thus in our iTEBD simulations, 
we never expect to see strict revivals. 

To illustrate these time scales in operation, we con¬ 
sider the occupation number, rii(x) = (x)ipi(x), for a 
fermionic mode on chain i, a simple measure of how the 
system departs from the initial state, for which rii{x) = 0. 
In Fig. 2 we present how rii(x) evolves with time for a 
quench to j± = 0.1. On the basis of our perturbative 
results for very small j± [30], we plot n(x) in units of 
j]_ for all four quenches presented. These four quenches 
correspond to four different chain lengths, R. 

We see that at short times, the results for rii{x)/j\ 
collapse onto a single curve as a function of t/t\. As 
time increases, the curves cease to track one another. 
The first to do this is the R = 4 curve, then the R = 
6 curve, and then finally the R = 8 curve. The time 
at which this happens corresponds, roughly, to tR, the 
scale on which the quench explores the finite length of 
the chain. We expect a small departure from this time 
scale because a finite j± will renormalize the quasiparticle 
velocity v = 1 in t R . We also see from the inset of Fig. 2 
that the evolution at longer times is no longer described 
by perturbation theory. 

In Fig. 3 we explore a quench to a j± which exceeds 
j c , the critical coupling for the 2 + 1 dimensional system. 
Such a quench is among the most challenging numerically 
as the population of higher energy chain states becomes 
significant. Concomitantly, the time evolution is most 
dependent on E c in this case. Ramped, rather than sud¬ 
den, quenches can be implemented with some possible 
advantages in this regard [37], though we have not yet 
explored this possibility. Nonetheless in Fig. 3 we see 
that for a given chain length, R, we can find cutoffs, E c 
such that the time evolution is converged. 

It is also possible to calculate postquench correla¬ 
tions between the chains. We show the nearest neigh¬ 
bor spin-spin correlation function as a function of time, 
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FIG. 4. Logarithm of return probability G(t), for R = 6 for 
j± = 0.1,0.5. Non-analytic behaviour is seen at short times 
for a quench to j± = 0.5. We find no non-analytic points 
for the corresponding quench to j± = 0.1, even at longer 
times up to t = £a = 10.0 (not shown). Inset: comparison 
of the infinite chain number system data with a system with 
N = 100 chains computed using TEBD. The first non-analytic 
point for the infinite cylinder forms the edge of a plateau, 
whereas for a finite number of chains it takes the form of a 
peak. 

(af(x,t)a? +1 (x,t)), for a selection of R and j± in the 
inset of Fig. 3. Our choice of j± > 0 favors antifer¬ 
romagnetic correlations, producing the overall negative 
sign. An expansion in small t shows that this quantity 
is proportional to j±t 2 allowing us to collapse the results 
onto a single curve at short times. Here we see signa¬ 
tures of both the tj x and Ir scales. In the inset we have 
marked the intrachain scattering time tj L , for the system 
with R = 8 and j± = 0.1. It is visible as the time that 
the j_ l = 0.1 and j± = 0.01 data begin to diverge. We 
also mark the time scale tn at which the data for chains 
with R = 8,j± = 0.01 begins to diverge from that of 
R = 10, Al = 0.01. 

To show that our method can handle non-trivial as¬ 
pects of quenching through the critical coupling of the 
coupled chain system, we search for non-analyticities 
in the Loschmidt echo as a function in time. The 
‘Loschmidt echo’ or overlap probability at a particular 
t is the modulus squared of the overlap between the ini¬ 
tial and time evolved state: 

G{t) = |(4' 0 |e- !;ff2D « lt |T 0 )| 2 (3) 

where \I/o is the ground state of the uncoupled chain sys¬ 
tem. In ID it is useful to define a per site rate func¬ 
tion, £(t) via G(t) = exp[— N£(t)]. Non-analyticities in 
£(t) have been interpreted as ‘dynamical phase transi¬ 
tions’, following an exact calculation of this quantity for 
the ID TFIM [38-40]. The general association of such 


non-analytic points with equilibrium critical phenomena 
is contested [41, 42], but we demonstrate analytically in 
low order perturbation theory[30] that for quenches to 
jj_ > 0.27 we expect non-analyticities in G(f). While 
this estimate for the value of j± is larger than j c - be¬ 
cause of the low order to which we took the computation 
- it does suggest that simple perturbation theory for the 
quantity G{t ) can be used to estimate the phase bound¬ 
aries in some 2D quantum systems. 

In Fig. 4 we plot logG(f) for a quench to j± = 0.5 - 
a value of j± where we should see non-analyticities. In 
2D this quantity scales with system volume RN, as does 
its ID counterpart [38]. It also scales with j' x . As ex¬ 
pected we find non-analytic behaviour for this quench, 
within the time window we are able to simulate, and see 
that the non-analyticity has the same qualitative struc¬ 
ture for both E c = 7 |A| and 8 |A|. For comparison we 
plot logG(f) for a quench to j± =0.1, where in contrast 
we find that this quantity is smooth within our simula¬ 
tion window. We remark that non-analyticities appear 
for the same quantity with j± = 0.2 (not plotted), just 
above j c = 0.185, but they first occur only at the edge of 
the attainable times with iTEBD. 

Finally we consider the case of finite length and open 
boundary conditions. The TEBD algorithm is slower by 
approximately a factor of N due to the loss of transla¬ 
tional invariance along the cylinder. We find negligible 
effect, for finite N > 10 and i away from the ends of 
the cylinder, on the results for local quantities such as 
rii(x) (up to the time scales we reach). However this is 
not true for the Loschmidt echo (a global measure), es¬ 
pecially when |jjJ > j c . The inset of Fig. 4 shows the 
difference between the iTEBD and N = 100 results for 
R = 6, jj l = 0.5. While there is excellent agreement up to 
t ~ 1a (not shown), afterwards there is a clear change in 
the non-analytic point structure. We also find that this 
effect is even more pronounced for very small R and large 
N (where our model reduces to a single ID TFIM), sug¬ 
gesting that boundary conditions have a non-negligible 
effect on the Loschmidt echo even for large systems. This 
last result has important consequences for possible exper¬ 
imental investigations. 

Conclusions: We have demonstrated a robust method 
to compute dynamical behaviour in 2D quantum 
(d=2+l) systems after a quench, which we intend to use 
to study other systems including coupled quantum wires 
(i.e. coupled Luttinger liquids) and Heisenberg chains. 
The algorithm should prove especially useful when inter¬ 
preting non-equilibrium cold atom [43, 44] and pump- 
probe experiments in the cuprates [45, 46]. 
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John Cardy, Fabian Essler, Andrew Goldsborough, Is¬ 
rael Klich, Anatoli Polkonikov, Rudolf Rorner and Steve 
Simons. This work was supported by the Engineering 
and Physical Sciences Research Council (grant number 
EP/L010623/1) and the US Department of Energy, Of- 
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AC02-98CH10886. 


SUPPLEMENTAL MATERIAL 


Eree Fermions 

In this section we describe our method applied to an 
exactly solvable quantum model in 2D. Consider a free 
Majorana field i/> = if >^, {4>{x, t), 4>(x', t)} = S(x-x’) with 
mass A, confined to a ring and with a Lorentz invariant 
action 


S = J dt J dxip^j^dfj, — A )ip, (4) 


where p = 0,1 refers to time and space coordinates, = 
•tpi'y 0 and 7 0 ’ 1 are suitable 2D Dirac matrices. The held 
has an expansion in fermionic modes [a n ,a^] = <5 m>n 

uj*al l e i(Un ~ xPn) ), 

^ = ■»+ 

wa t e i{fc»-xrd) i ( 5 ) 


Here w = e*’ 1 '/ 4 , and e„ = Acoshd„, p n = Asinh0„. 
The momentum can take discrete values p n = 2irn/R for 
integer n, and e n = a/A 2 +p 2 . Note that in this case 
we are not mapping from a putative spin system using 
the Jordan-Wigner transformation, so there is no sepa¬ 
ration of the spectrum into Ramond and Neveu-Sclrwarz 
sectors. Using the mode expansion we obtain the chain 
Hamiltonian 


Hid/ = ^ U 




( 6 ) 


with chain index £. We can build a 2D quantum system 
from an array of N of these chains coupled with a nearest 
neighbor interaction (here we assume our 2D system is a 
torus) 


Hint = 


\ - A , ! 

— ~\ a n,e a n ,£+1 

e,n n 


+ H.C.) 


( 7 ) 


index n along the chains). 


_ 


-^E- 


j /I ^k m ,k m r 


ik m C J 

a i ; 


2 At i 


Hf r ee — ^ COS km ] ah m Ur* 






( 8 ) 

(9) 

( 10 ) 


The Hamiltonian is diagonal in n and in : if so desired it 
can be treated as either a set of N ID uncoupled bands 
indexed by m, or infinitely many ID bands, indexed by n. 
In either case the system can then be treated by standard 
ID MPS methods, because the lack of coupling between 
bands means that one needs to keep only those eigen¬ 
states from the spectrum of Hyd, e that are present in the 
initial state of the chain array. Beyond this the cutoff 
E c does not play a role. For example consider the initial 
state 


l$> n ^( |0)2i 

2=0 


l n = 0 > 2 i+l + l n = 0 ) 2 il°)i+l)’ 


( 11 ) 


in which alternating pairs of chains are entangled, with a 
superposition of ground states (( 0 )^ and lowest excited 
states (|n = 0)^) on the chains. Evolving this state under 
Hf ree does not involve any other states from the chain 
spectrum. For such an evolution the return probability 
can be calculated exactly using a determinant method, 
yielding 


G(t) = |det(M + (1 — M)Hexp{ — iht}R) 
M = diag(0,1,0,1,0, •••), 

/l 1 \ 


( 12 ) 


R = 


1 -1 


1 1 
1 -1 


V 


■J 


( A \ 


h = 


-t± A 


(13) 


with N x N matrices. We compare this result (evaluated 
for a torus of N = 800 chains) with the iTEBD result 
computed with our code implementing Hf ree in Fig. 5. 
Note that in this case a change in A can be absorbed into 
a simultaneous rescaling of t± and t, so the hopping t± 
sets the time scale. 


where the coupling strength is t±. This chain array sys¬ 
tem Hfree = Yle Hiv,e + H lnt can be solved trivially by 
Fourier transformation from chain index £ to momentum 
k m = 2 tt m/N (note this is transverse to the momentum 


Further details of method 

The structure of the spectrum of the continuum limit 
Ising chain is detailed in Ref. [36] but we give a brief de- 
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FIG. 5. Logarithm of return probability, G(t), for the free 
fermion quench described in the text, with hopping parameter 
t± = 0.5. Both the exact result, using a determinant method 
with N = 800 chains, and the iTEBD result are plotted. In¬ 
set: enlarged region showing the determinant calculation and 
iTEBD with different s m in parameter. 

scription here for convenience. The spectrum splits into 
two sectors, Neveu-Schwarz (NS) and Rarnond (RM). 
The energy of a state with a particular fermion configu¬ 
ration is given by 


E({n s }) =E S +J2 E n s (14) 

{n fl } 



where s = NV or RM and E s is the vacuum energy (dif¬ 
ferent in the two sectors). For the disordered phase of 
a chain, A < 0, states with even numbers of particles 
(including the 0 particle vacuum state) are in the NS 
sector (n s £ Z), while odd particle states are in the Ra- 
mond sector (n s £ Z/2) The spin operator, <r z , is off 
diagonal in sector, so that on an individual chain RM 
states are only scattered into NS states and vice-versa by 
the J_l term. This fact that makes perturbative calcu¬ 
lations significantly easier. As a consequence the overall 
sector (whether there is a odd or even number of Ra- 
mond chains) is conserved by H 2 dqi- The total sector 
and momentum (along the chain direction) for two chains 
is also conserved by the two chain time evolution oper¬ 
ator. These conservation laws are useful for performing 
matrix operations by sub blocks. 

The fermionic representation is symmetric with respect 
to the spin direction, as is our pre-quench state, and 
i^ 2 DQi itself does not break this symmetry. Hence the 
local magnetization (crf(x)) is always zero. 

We implement iTEBD [31] using the alteration due 
to Hastings [47] that improves numerical robustness by 



FIG. 6. Convergence in s m w for R = 10, J± = 0.1, using 
a smallest retained singular value criterion, s m in (the largest 
bond dimension used during the calculation is given as Xf- 
Also shown for comparison is data collected using a fixed bond 
dimension, %■ 

removing the need to divide by very small singular values. 
Before computing expectation values the iTEBD transfer 
matrix must be ‘orthogonalised’ as described in Refs. [14, 
32, 48]. 

Instead of imposing a fixed bond dimension, we per¬ 
form TEBD and iTEBD using a cutoff on the minimum 
singular value that is retained, s m ; n . This translates to 
a minimum eigenvalue of the reduced density matrix, 
Pmin = Smin- The advantage over fixed bond dimen¬ 
sion is that when working with large matrices, one may 
drop unwanted singular values as individual sub blocks 
are processed, as opposed to recording the full results of 
the singular value decomposition and sorting by magni¬ 
tude before truncating. In principle the maximum value 
of Se that can be captured in this way is —2 log s m in- In 
Fig. 6 we give an example of convergence in this param¬ 
eter and compare with a version of the algorithm that 
uses a fixed bond dimension throughout. 

We mainly use the second order Trotter decomposition 
in this work, and find it sufficient for the time scales we 
are able to reach before bond dimension, y, becomes the 
limiting factor. For very weak quenches however, a first 
order Trotter decomposition still gives good results, even 
to quite large times. This is because the leading error 
term in the decomposition is proportional to J±. A more 
sensitive test is whether the sharp non-analytic features 
in the return probability (Loschmidt echo) for j± = 0.5 
are well converged in the time step r. Figure 7 shows 
the logarithm of the return probability for r = 0.0081 a 
and 0.002£a using second order Trotter decompositions. 
There is negligible difference except in the immediate ap¬ 
proach to the non-analytic point, and even in this case 
the difference is less than 1%. 
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FIG. 7. Convergence in the time step r (in units of 1 a). Left 
panel: In the perturbative (j± -C 1) limit results converge 
extremely quickly as a function of t, for both 1st order (Tl) 
and 2nd order (T 2) Trotter decompositions. Right panel: At 
the other extreme, for j± > j c we see that the sensitive place¬ 
ment of non-analytic points is converged for sufficiently small 
r. 


Discussion of the effects of the cutoff 


TABLE I. The number of chain states kept (including the 
ground state) for various combinations of R and E c (A = — 1). 


R 

Ec/\A\ 

4 5 6 8 10 

2 

5 5 5 12 19 

4 

11 14 19 43 90 

6 

19 33 52 124 

8 

29 55 103 308 

10 

44 97 172 

12 

63 


|A|i? <C 1 only the chain ground state and lowest excited 
state survive with the energies of higher excited states 
scaling as multiples of (|A|i?) -1 . For A < 0 the ground 
and first excited states are the Neveu-Schwarz vacuum 
and the zero momentum single particle Rarnond states 
respectively. In this case the Hamiltonian of the coupled 
chain system becomes 


lim H 2 dqi = T 
R-> o 


( A + E r m 

0 ^ 

K 0 

Ejsss ) 


-\-J±R 


0 M 
M 0 


0 M 
M 0 


i +1 

(16) 


where M = (RM, k = 0| a z |NS,vac) € R. This can be 
written as the Hamiltonian of a single ID lattice quantum 
Ising chain (up to some unimportant constants): 

hm H 2 bqi = L/idqi = 5Z °f+ 1 ) > ( 17 ) 

i 


where a x ,a z are the usual Pauli matrices and we make 
the identifications 


h= l A l + ^M-^NS ; (18) 

J = J ± RM 2 . (19) 


This is a useful check of the code, as it is easy establish if 
one recovers the correct ID behaviour, including the 1 + 1 
dimensional phase transition when h = J. For example 
at R = 1, and by using E c to restrict the number of chain 
states to two, we are able to successfully reproduce the 
predicted positions of non-analytic points, for quenches 
of the ID quantum Ising model through its critical point 

[38]. 


Perturbation Theory 


The number of chain states increases approximately 
exponentially with E c and R for E c > A. In Table I we 
show the number of chain states kept for a range of chain 
lengths and cutoffs. It is clear that changing E c from 
4 | A | to 8 | A | (for example) has a much more dramatic 
effect on the number of included states at R = 8 than 
at R = 4. However with quenches below the critical 
coupling J c , we find that in general we see very little 
difference between E c = 6 |A| and ^ 8 |A| for R > 4. 


Small R limit 

For a single chain there should be a crossover to effec¬ 
tive OD (0+1 dimensional) behaviour when the correla¬ 
tion length is of order the system size, | A |~ 1 ~ R. When 


In the limit of small interchain coupling \J±/J C \ 1 

a perturbative expansion is appropriate. We use uni¬ 
tary perturbation theory, following Ref. [49], in order to 
avoid spurious secular terms that grow in time without 
bound. The expectation of an operator A at time t after 
the quench (assuming that the pre quench Hamiltonian 
commutes with it, [Hid,*, A] = 0) is given to order J\ by 


(A(t)) = (4(0)) + 4Ji /°°da;FH Sln ~y /2) 

J — OO ^ 

+( W )=^|(4/ 0 |^ ( 7^f +1 |ct>)| 2 

<l> i 

x ($| A |$><S(w-(£*-£*„)), 


( 20 ) 

( 21 ) 

( 22 ) 


where T) is a tensor product of unperturbed chain states, 
|<^i) ® |^ 2 ) ( 8 > ■ ■ ■ < 8 > I^jv)- In this work we choose the pre 
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FIG. 8. Fermion occupation number, rii(x), at position x on 
chain i, scaled by the interchain coupling, J± = 0.01, squared. 
Differences as a function of cutoff are negligible for E c > 4 | A|. 
Dashed line: the perturbative result for R = 10. 


quench state to be the ground state of the uncoupled 
chain system, 

i*o>=ni°‘>> ( 2s ) 


where |0*) is the ground state of chain i. For the operator 
A we consider the occupation number for a fermion 
on chain i with momentum k (along the chain). With 
these choices the only states that contribute to F(w) for 
the quantum Ising system are tensor products of chain 
vacuum states, with a single nearest neighbor pair of 
chain excited states: 


1$+) = • • • 0 |0j_i) 0 10j) 0 |0 i+ i) 0 |0j+ 2 ) 0 • • • , 

|$_) = • • ■ 0 |0j_ 2 ) 0 |0i_i) 0 |</>i) 0 |0,+i) 0 • • ■ . (24) 


Using these states we find 

f 

F(uj)=R 2 ^2 <0*1 <7? 10,) (0i +t \<r* i+ s\4>i + ,) 


q>i+s,q>i 
s=± 1 


x 6(u - (Efc + E </>i+s - 2 E 0 )) ( 4>i \nqfe | fa) 


(25) 


where E^ is the energy of the (unperturbed) chain state 
| <f>). The restriction on the sum indicates that the mo¬ 
menta of the chain states | <f>i ), must sum to zero. 

For the expectation at time t we obtain 


/ 

K,fc(i)) = {2J ± R) 2 Y^ | (OiK |^»)(0 i+s | of +s | (j> i+a ) 

S =± 1 


X 


sin 2 (tF; i+Sii /2) 

Ef +S ^ 


{4>i I l R'i,k \4*i) • 


(26) 


Here E l+ai = E ^ + E ( j >i+s — 2 Eq and the restriction on 
the sum is as above. The sum over s will contribute a 
simple factor of 2 unless the system is an open cylinder 
and i = 1, N, in which case one of the sums vanishes due 
to the missing nearest neighbor. 

We now make some remarks about the result at or¬ 
der Jj_. For all i. k Eq. 26 is a sum of oscillatory terms 
with no ‘decay’ even in the thermodynamic limit. There 
are no boundary effects, excepting the trivial factor of 2 
described above due to the different number of nearest 
neighbors. With disordered chains, for which the chain 
ground state is in the Neveu-Schwarz sector, only excited 
chain states of the Ramond sector will contribute. Con¬ 
sequently at this order (nq^) will be zero for half integer 
momenta, k = 2i r(n + l/2)/i?. 

Using the above we calculate the occupation numbers 
perturbatively for a large range of k by evaluating the 
sums over states numerically. From these results it is 
simple to to find the position space occupation: 

f R 

Rrii(x) = / dx 7ij{x) = (27) 

Jo 


using translational invariance along the chains, provided 
the rii } k drop off sufficiently rapidly with k. We show 
our results for J_l = 0.01 together with the perturbative 
curve (dashed curve) for R = 10 in Fig. 8. The agree¬ 
ment is excellent at short times and still very good at 
longer times for rii(x). 


Estimate of Critical Coupling, J c 

In this portion of the supplemental material we ar¬ 
gue that the appearance of non-analyticities in the return 
probability post-quench can be used to estimate the value 
of the critical coupling J c marking the phase transition 
(although see discussion below). To this end we employ 
the observation of [38] that in quenching from an initial 
value of the coupling, J± % . to a final value of the cou¬ 
pling, Jl/, non-analyticities appear in the return prob¬ 
ability whenever a fermionic mode (of the post-quench 
Hamiltonian) has an occupation of at least 1/2, i.e. ap¬ 
pears with an occupation corresponding to either infinite 
or negative temperatures. In the case considered in [38], 
these athermal occupations only occurred if the coupling 
crossed a phase boundary. The appearance of the non- 
analyticities can then be used to estimate the location of 
these boundaries. However it is at least possible that in 
general interacting models (such as the XXZ spin chain 
considered in [41, 42], such occupations can be induced 
without crossing a phase boundary. We will use low order 
perturbation theory to estimate the coupling J c at which 
athermal mode occupations appear, noting that because 
of the work of [41, 42] that this J c may not correspond 
to the critical coupling determining the equilibrium phase 
transition in the two dimensional quantum Ising model. 
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The modes that we will consider in this argument take 
the form 

4 + " ' ( 28 ) 

where A 4 fe is an operator on the j-th chain that creates 
a fermion with momentum k x along the chain. We then 
want to find the minimum value of J± such that 

(nk x ,k y ) = (i\4>i„k y ( J -L)' l Pk :c ,k y ( J ±)\i) = 1 / 2 , 

where |i) is the initial state of the quench (here the 
ground state of the system for J± = 0). The mode 
for which this will first occur is (k x ,k y ) = (fcx.mm, 0) = 
(2-7T /R, 0) as this is the mode with the lowest energy that 
couples to the perturbation, and so is easiest to drive 
athermal. 

To compute mi „,o) we expand it in terms of the 
eigenstates {| s)} of the post-quench Hamiltonian 

( n k x , min ,o) = l< s N>| 2 ( s K*, m in,o|s>- (29) 


We will suppose that this sum is dominated by states 
involving at most one fermion on any given chain. The 
only such state |s) that then contributes to this sum is 

| k x ,min , 9, k XiTn i n , 

While we cannot write down an exact expression for this 
state as a function of Jj_, we are able to write down to 
second order the contribution to this state coming from 
the Jj_ = 0 vacuum |0) - the only part that matters in 
computing the overlap (s|i). To second order we have 

\k x ,mini 0; —kx^min, 0) = mirt ,o(®)' i /-fc x , m i„ ,o(®) 1^) 

+ -^2 (* - 1^2 - b C , ('li))|0) ■+ • (30) 

We use here the conventions and notation of [36]. In 
particular a = 1.35783834 ... and = A 2 + k x . Thus 
to second order in J j_ we have 




|A| 1/2 d 4 Ji 4|A| 1//2 ct 4 

n, x, min ry, x ,mzn 


(31) 


For R —> oo the value of j± = J_l/|A|'^ 4 at which 

( n fc a ,,™iT,,o) = 1/2 is j _l = 0.27119_ While this value 

of j _l is considerably larger than the value of j c = 0.185 
given in [21] for where the equilibrium phase transition 
occurs, we can at least partially attribute this difference 
to the neglect both of terms higher order in J± in this 
computation as well as Jj_ = 0 states involving more than 
one fermion per chain. 


* andrew.james@ucl.ac.uk 
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